home *** CD-ROM | disk | FTP | other *** search
/ EuroCD 3 / EuroCD 3.iso / Programming / Python1.4_Source / Objects / complexobject.c < prev    next >
C/C++ Source or Header  |  1998-06-24  |  12KB  |  628 lines
  1. /* Complex object implementation */
  2.  
  3. /* Borrows heavily from floatobject.c */
  4.  
  5. #ifndef WITHOUT_COMPLEX
  6.  
  7. #include "allobjects.h"
  8. #include "modsupport.h"
  9. #include "protos/complexobject_protos.h"
  10. #include <errno.h>
  11. #include "mymath.h"
  12.  
  13. #ifdef i860
  14. /* Cray APP has bogus definition of HUGE_VAL in <math.h> */
  15. #undef HUGE_VAL
  16. #endif
  17.  
  18. #ifdef HUGE_VAL
  19. #define CHECK(x) if (errno != 0) ; \
  20.     else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
  21.     else errno = ERANGE
  22. #else
  23. #define CHECK(x) /* Don't know how to check */
  24. #endif
  25.  
  26. #ifdef HAVE_LIMITS_H
  27. #include <limits.h>
  28. #endif
  29.  
  30. #ifndef LONG_MAX
  31. #define LONG_MAX 0X7FFFFFFFL
  32. #endif
  33.  
  34. #ifndef LONG_MIN
  35. #define LONG_MIN (-LONG_MAX-1)
  36. #endif
  37.  
  38. #ifdef __NeXT__
  39. #ifdef __sparc__
  40. /*
  41.  * This works around a bug in the NS/Sparc 3.3 pre-release
  42.  * limits.h header file.
  43.  * 10-Feb-1995 bwarsaw@cnri.reston.va.us
  44.  */
  45. #undef LONG_MIN
  46. #define LONG_MIN (-LONG_MAX-1)
  47. #endif
  48. #endif
  49.  
  50. #if !defined(__STDC__) && !defined(macintosh)
  51. extern double fmod PROTO((double, double));
  52. extern double pow PROTO((double, double));
  53. #endif
  54.  
  55.  
  56. /* elementary operations on complex numbers */
  57.  
  58. static int c_error;
  59. static Py_complex c_1 = {1., 0.};
  60.  
  61. Py_complex c_sum(a,b)
  62.     Py_complex a,b;
  63. {
  64.     Py_complex r;
  65.     r.real = a.real + b.real;
  66.     r.imag = a.imag + b.imag;
  67.     return r;
  68. }
  69.  
  70. Py_complex c_diff(a,b)
  71.     Py_complex a,b;
  72. {
  73.     Py_complex r;
  74.     r.real = a.real - b.real;
  75.     r.imag = a.imag - b.imag;
  76.     return r;
  77. }
  78.  
  79. Py_complex c_neg(a)
  80.     Py_complex a;
  81. {
  82.     Py_complex r;
  83.     r.real = -a.real;
  84.     r.imag = -a.imag;
  85.     return r;
  86. }
  87.  
  88. Py_complex c_prod(a,b)
  89.     Py_complex a,b;
  90. {
  91.     Py_complex r;
  92.     r.real = a.real*b.real - a.imag*b.imag;
  93.     r.imag = a.real*b.imag + a.imag*b.real;
  94.     return r;
  95. }
  96.  
  97. Py_complex c_quot(a,b)
  98.     Py_complex a,b;
  99. {
  100.     Py_complex r;
  101.     double d = b.real*b.real + b.imag*b.imag;
  102.     if (d == 0.)
  103.         c_error = 1;
  104.     r.real = (a.real*b.real + a.imag*b.imag)/d;
  105.     r.imag = (a.imag*b.real - a.real*b.imag)/d;
  106.     return r;
  107. }
  108.  
  109. Py_complex c_pow(a,b)
  110.     Py_complex a,b;
  111. {
  112.     Py_complex r;
  113.     double vabs,len,at,phase;
  114.     if (b.real == 0. && b.imag == 0.) {
  115.         r.real = 1.;
  116.         r.imag = 0.;
  117.     }
  118.     else if (a.real == 0. && a.imag == 0.) {
  119.         if (b.imag != 0. || b.real < 0.)
  120.             c_error = 2;
  121.         r.real = 0.;
  122.         r.imag = 0.;
  123.     }
  124.     else {
  125.         vabs = hypot(a.real,a.imag);
  126.         len = pow(vabs,b.real);
  127.         at = atan2(a.imag, a.real);
  128.         phase = at*b.real;
  129.         if (b.imag != 0.0) {
  130.             len /= exp(at*b.imag);
  131.             phase += b.imag*log(vabs);
  132.         }
  133.         r.real = len*cos(phase);
  134.         r.imag = len*sin(phase);
  135.     }
  136.     return r;
  137. }
  138.  
  139. static Py_complex c_powu(x, n)
  140.     Py_complex x;
  141.     long n;
  142. {
  143.     Py_complex r, p;
  144.     long mask = 1;
  145.     r = c_1;
  146.     p = x;
  147.     while (mask > 0 && n >= mask) {
  148.         if (n & mask)
  149.             r = c_prod(r,p);
  150.         mask <<= 1;
  151.         p = c_prod(p,p);
  152.     }
  153.     return r;
  154. }
  155.  
  156. static Py_complex c_powi(x, n)
  157.     Py_complex x;
  158.     long n;
  159. {
  160.     Py_complex cn;
  161.  
  162.     if (n > 100 || n < -100) {
  163.         cn.real = (double) n;
  164.         cn.imag = 0.;
  165.         return c_pow(x,cn);
  166.     }
  167.     else if (n > 0)
  168.         return c_powu(x,n);
  169.     else
  170.         return c_quot(c_1,c_powu(x,-n));
  171.  
  172. }
  173.  
  174. PyObject *
  175. PyComplex_FromCComplex(cval)
  176.     Py_complex cval;
  177. {
  178.     register complexobject *op =
  179.         (complexobject *) malloc(sizeof(complexobject));
  180.     if (op == NULL)
  181.         return err_nomem();
  182.     op->ob_type = &Complextype;
  183.     op->cval = cval;
  184.     NEWREF(op);
  185.     return (object *) op;
  186. }
  187.  
  188. PyObject *
  189. PyComplex_FromDoubles(real, imag)
  190.     double real, imag;
  191. {
  192.     Py_complex c;
  193.     c.real = real;
  194.     c.imag = imag;
  195.     return PyComplex_FromCComplex(c);
  196. }
  197.  
  198. double
  199. PyComplex_RealAsDouble(op)
  200.     PyObject *op;
  201. {
  202.     if (PyComplex_Check(op)) {
  203.         return ((PyComplexObject *)op)->cval.real;
  204.     } else {
  205.         return PyFloat_AsDouble(op);
  206.     }
  207. }
  208.  
  209. double
  210. PyComplex_ImagAsDouble(op)
  211.     PyObject *op;
  212. {
  213.     if (PyComplex_Check(op)) {
  214.         return ((PyComplexObject *)op)->cval.imag;
  215.     } else {
  216.         return 0.0;
  217.     }
  218. }
  219.  
  220. Py_complex
  221. PyComplex_AsCComplex(op)
  222.     PyObject *op;
  223. {
  224.     Py_complex cv;
  225.     if (PyComplex_Check(op)) {
  226.         return ((PyComplexObject *)op)->cval;
  227.     } else {
  228.         cv.real = PyFloat_AsDouble(op);
  229.         cv.imag = 0.;
  230.         return cv;
  231.     }   
  232. }
  233.  
  234. static void
  235. complex_dealloc(op)
  236.     object *op;
  237. {
  238.     DEL(op);
  239. }
  240.  
  241.  
  242. static void
  243. complex_buf_repr(buf, v)
  244.     char *buf;
  245.     complexobject *v;
  246. {
  247.     if (v->cval.real == 0.)
  248.         sprintf(buf, "%.12gj", v->cval.imag);
  249.     else
  250.         sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
  251. }
  252.  
  253. static int
  254. complex_print(v, fp, flags)
  255.     complexobject *v;
  256.     FILE *fp;
  257.     int flags; /* Not used but required by interface */
  258. {
  259.     char buf[100];
  260.     complex_buf_repr(buf, v);
  261.     fputs(buf, fp);
  262.     return 0;
  263. }
  264.  
  265. static object *
  266. complex_repr(v)
  267.     complexobject *v;
  268. {
  269.     char buf[100];
  270.     complex_buf_repr(buf, v);
  271.     return newstringobject(buf);
  272. }
  273.  
  274. static int
  275. complex_compare(v, w)
  276.     complexobject *v, *w;
  277. {
  278. /* Note: "greater" and "smaller" have no meaning for complex numbers,
  279.    but Python requires that they be defined nevertheless. */
  280.     Py_complex i, j;
  281.     i = v->cval;
  282.     j = w->cval;
  283.     if (i.real == j.real && i.imag == j.imag)
  284.        return 0;
  285.     else if (i.real != j.real)
  286.        return (i.real < j.real) ? -1 : 1;
  287.     else
  288.        return (i.imag < j.imag) ? -1 : 1;
  289. }
  290.  
  291. static long
  292. complex_hash(v)
  293.     complexobject *v;
  294. {
  295.     double intpart, fractpart;
  296.     int expo;
  297.     long x;
  298.     /* This is designed so that Python numbers with the same
  299.        value hash to the same value, otherwise comparisons
  300.        of mapping keys will turn out weird */
  301.  
  302. #ifdef MPW /* MPW C modf expects pointer to extended as second argument */
  303. {
  304.     extended e;
  305.     fractpart = modf(v->cval.real, &e);
  306.     intpart = e;
  307. }
  308. #else
  309.     fractpart = modf(v->cval.real, &intpart);
  310. #endif
  311.  
  312.     if (fractpart == 0.0) {
  313.         if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
  314.             /* Convert to long int and use its hash... */
  315.             object *w = dnewlongobject(v->cval.real);
  316.             if (w == NULL)
  317.                 return -1;
  318.             x = hashobject(w);
  319.             DECREF(w);
  320.             return x;
  321.         }
  322.         x = (long)intpart;
  323.     }
  324.     else {
  325.         fractpart = frexp(fractpart, &expo);
  326.         fractpart = fractpart*2147483648.0; /* 2**31 */
  327.         x = (long) (intpart + fractpart) ^ expo; /* Rather arbitrary */
  328.     }
  329.     if (x == -1)
  330.         x = -2;
  331.     return x;
  332. }
  333.  
  334. static object *
  335. complex_add(v, w)
  336.     complexobject *v;
  337.     complexobject *w;
  338. {
  339.     return newcomplexobject(c_sum(v->cval,w->cval));
  340. }
  341.  
  342. static object *
  343. complex_sub(v, w)
  344.     complexobject *v;
  345.     complexobject *w;
  346. {
  347.     return newcomplexobject(c_diff(v->cval,w->cval));
  348. }
  349.  
  350. static object *
  351. complex_mul(v, w)
  352.     complexobject *v;
  353.     complexobject *w;
  354. {
  355.     return newcomplexobject(c_prod(v->cval,w->cval));
  356. }
  357.  
  358. static object *
  359. complex_div(v, w)
  360.     complexobject *v;
  361.     complexobject *w;
  362. {
  363.     Py_complex quot;
  364.     c_error = 0;
  365.     quot = c_quot(v->cval,w->cval);
  366.     if (c_error == 1) {
  367.         err_setstr(ZeroDivisionError, "complex division");
  368.         return NULL;
  369.     }
  370.     return newcomplexobject(quot);
  371. }
  372.  
  373. static object *
  374. complex_remainder(v, w)
  375.     complexobject *v;
  376.     complexobject *w;
  377. {
  378.         Py_complex div, mod;
  379.     c_error = 0;
  380.     div = c_quot(v->cval,w->cval); /* The raw divisor value. */
  381.     if (c_error == 1) {
  382.         err_setstr(ZeroDivisionError, "complex remainder");
  383.         return NULL;
  384.     }
  385.     div.real = floor(div.real); /* Use the floor of the real part. */
  386.     div.imag = 0.0;
  387.     mod = c_diff(v->cval, c_prod(w->cval, div));
  388.  
  389.     return newcomplexobject(mod);
  390. }
  391.  
  392.  
  393. static object *
  394. complex_divmod(v, w)
  395.     complexobject *v;
  396.     complexobject *w;
  397. {
  398.         Py_complex div, mod;
  399.     PyObject *d, *m, *z;
  400.     c_error = 0;
  401.     div = c_quot(v->cval,w->cval); /* The raw divisor value. */
  402.     if (c_error == 1) {
  403.         err_setstr(ZeroDivisionError, "complex divmod()");
  404.         return NULL;
  405.     }
  406.     div.real = floor(div.real); /* Use the floor of the real part. */
  407.     div.imag = 0.0;
  408.     mod = c_diff(v->cval, c_prod(w->cval, div));
  409.     d = newcomplexobject(div);
  410.     m = newcomplexobject(mod);
  411.     z = mkvalue("(OO)", d, m);
  412.     Py_XDECREF(d);
  413.     Py_XDECREF(m);
  414.     return z;
  415. }
  416.  
  417. static object *
  418. complex_pow(v, w, z)
  419.     complexobject *v;
  420.     object *w;
  421.     complexobject *z;
  422. {
  423.     Py_complex p;
  424.     Py_complex exponent;
  425.     long int_exponent;
  426.  
  427.      if ((object *)z!=None) {
  428.         err_setstr(ValueError, "complex modulo");
  429.         return NULL;
  430.     }
  431.  
  432.     c_error = 0;
  433.     exponent = ((complexobject*)w)->cval;
  434.     int_exponent = (long)exponent.real;
  435.     if (exponent.imag == 0. && exponent.real == int_exponent)
  436.         p = c_powi(v->cval,int_exponent);
  437.     else
  438.         p = c_pow(v->cval,exponent);
  439.  
  440.     if (c_error == 2) {
  441.         err_setstr(ValueError, "0.0 to a negative or complex power");
  442.         return NULL;
  443.     }
  444.  
  445.     return newcomplexobject(p);
  446. }
  447.  
  448. static object *
  449. complex_neg(v)
  450.     complexobject *v;
  451. {
  452.     Py_complex neg;
  453.     neg.real = -v->cval.real;
  454.     neg.imag = -v->cval.imag;
  455.     return newcomplexobject(neg);
  456. }
  457.  
  458. static object *
  459. complex_pos(v)
  460.     complexobject *v;
  461. {
  462.     INCREF(v);
  463.     return (object *)v;
  464. }
  465.  
  466. static object *
  467. complex_abs(v)
  468.     complexobject *v;
  469. {
  470.     return newfloatobject(hypot(v->cval.real,v->cval.imag));
  471. }
  472.  
  473. static int
  474. complex_nonzero(v)
  475.     complexobject *v;
  476. {
  477.     return v->cval.real != 0.0 && v->cval.imag != 0.0;
  478. }
  479.  
  480. static int
  481. complex_coerce(pv, pw)
  482.     object **pv;
  483.     object **pw;
  484. {
  485.     Py_complex cval;
  486.     cval.imag = 0.;
  487.     if (is_intobject(*pw)) {
  488.         cval.real = (double)getintvalue(*pw);
  489.         *pw = newcomplexobject(cval);
  490.         INCREF(*pv);
  491.         return 0;
  492.     }
  493.     else if (is_longobject(*pw)) {
  494.         cval.real = dgetlongvalue(*pw);
  495.         *pw = newcomplexobject(cval);
  496.         INCREF(*pv);
  497.         return 0;
  498.     }
  499.     else if (is_floatobject(*pw)) {
  500.         cval.real = getfloatvalue(*pw);
  501.         *pw = newcomplexobject(cval);
  502.         INCREF(*pv);
  503.         return 0;
  504.     }
  505.     return 1; /* Can't do it */
  506. }
  507.  
  508. static object *
  509. complex_int(v)
  510.     object *v;
  511. {
  512.     err_setstr(TypeError,
  513.            "can't convert complex to int; use e.g. int(abs(z))");
  514.     return NULL;
  515. }
  516.  
  517. static object *
  518. complex_long(v)
  519.     object *v;
  520. {
  521.     err_setstr(TypeError,
  522.            "can't convert complex to long; use e.g. long(abs(z))");
  523.     return NULL;
  524. }
  525.  
  526. static object *
  527. complex_float(v)
  528.     object *v;
  529. {
  530.     err_setstr(TypeError,
  531.            "can't convert complex to float; use e.g. abs(z)");
  532.     return NULL;
  533. }
  534.  
  535.  
  536. static object *
  537. complex_new(self, args)
  538.     object *self;
  539.     object *args;
  540. {
  541.     Py_complex cval;
  542.  
  543.     cval.imag = 0.;
  544.     if (!PyArg_ParseTuple(args, "d|d", &cval.real, &cval.imag))
  545.         return NULL;
  546.     return newcomplexobject(cval);
  547. }
  548.  
  549. static object *
  550. complex_conjugate(self)
  551.     object *self;
  552. {
  553.     Py_complex c;
  554.     c = ((complexobject *)self)->cval;
  555.     c.imag = -c.imag;
  556.     return newcomplexobject(c);
  557. }
  558.  
  559. static PyMethodDef complex_methods[] = {
  560.     {"conjugate",    (PyCFunction)complex_conjugate,    1},
  561.     {NULL,        NULL}        /* sentinel */
  562. };
  563.  
  564.  
  565. static object *
  566. complex_getattr(self, name)
  567.     complexobject *self;
  568.     char *name;
  569. {
  570.     Py_complex cval;
  571.     if (strcmp(name, "real") == 0)
  572.         return (object *)newfloatobject(self->cval.real);
  573.     else if (strcmp(name, "imag") == 0)
  574.         return (object *)newfloatobject(self->cval.imag);
  575.     else if (strcmp(name, "conj") == 0) {
  576.         cval.real = self->cval.real;
  577.         cval.imag = -self->cval.imag;
  578.         return (object *)newcomplexobject(cval);
  579.     }
  580.     return findmethod(complex_methods, (object *)self, name);
  581. }
  582.  
  583. static number_methods complex_as_number = {
  584.     (binaryfunc)complex_add, /*nb_add*/
  585.     (binaryfunc)complex_sub, /*nb_subtract*/
  586.     (binaryfunc)complex_mul, /*nb_multiply*/
  587.     (binaryfunc)complex_div, /*nb_divide*/
  588.     (binaryfunc)complex_remainder,    /*nb_remainder*/
  589.     (binaryfunc)complex_divmod,    /*nb_divmod*/
  590.     (ternaryfunc)complex_pow, /*nb_power*/
  591.     (unaryfunc)complex_neg, /*nb_negative*/
  592.     (unaryfunc)complex_pos, /*nb_positive*/
  593.     (unaryfunc)complex_abs, /*nb_absolute*/
  594.     (inquiry)complex_nonzero, /*nb_nonzero*/
  595.     0,        /*nb_invert*/
  596.     0,        /*nb_lshift*/
  597.     0,        /*nb_rshift*/
  598.     0,        /*nb_and*/
  599.     0,        /*nb_xor*/
  600.     0,        /*nb_or*/
  601.     (coercion)complex_coerce, /*nb_coerce*/
  602.     (unaryfunc)complex_int, /*nb_int*/
  603.     (unaryfunc)complex_long, /*nb_long*/
  604.     (unaryfunc)complex_float, /*nb_float*/
  605.     0,        /*nb_oct*/
  606.     0,        /*nb_hex*/
  607. };
  608.  
  609. typeobject Complextype = {
  610.     OB_HEAD_INIT(&Typetype)
  611.     0,
  612.     "complex",
  613.     sizeof(complexobject),
  614.     0,
  615.     (destructor)complex_dealloc,    /*tp_dealloc*/
  616.     (printfunc)complex_print,    /*tp_print*/
  617.     (getattrfunc)complex_getattr,    /*tp_getattr*/
  618.     0,                /*tp_setattr*/
  619.     (cmpfunc)complex_compare,    /*tp_compare*/
  620.     (reprfunc)complex_repr,        /*tp_repr*/
  621.     &complex_as_number,            /*tp_as_number*/
  622.     0,                /*tp_as_sequence*/
  623.     0,                /*tp_as_mapping*/
  624.     (hashfunc)complex_hash,     /*tp_hash*/
  625. };
  626.  
  627. #endif
  628.